%% fig s4e single neuron place field correlation directly
load( '/media/Big1/15tracks/HeadDirection/head_direction_selectivity_byPlfieldCorr.mat');



ccpara=[];ccorth=[];p1=[];
for irat=1:5
    for idir=1:2
        for ic=1:length(res(irat,idir).res1)
            c1=mean(res(irat,idir).res1(ic).datacc{1}{4});
            c2=mean(res(irat,idir).res1(ic).datacc{1}{3});
%             c1=(res(irat,idir).res1(ic).datacc{1}{4});
%             c2=(res(irat,idir).res1(ic).datacc{1}{3});
            ccpara=[ccpara;c1];
            ccorth=[ccorth;c2];
%                p1=[p1;ranksum(c1,c2)];
        end
    end
end


[m,se]=mean_se_cell({ccpara,ccorth});
pval=ranksum(ccpara,ccorth);
figure;
% errorbarformal(1:2,m,se,{'parallel','orthogonal'});
[y,group]=cell2group({ccpara,ccorth});
plot_errorbar_with_data(y,group,{'parallel','orthogonal'});
sigstar([1,2],pval);
ylabel('Correlation');
title('single neuron place field');
setaxisformal

% savepdf(gcf,'/media/Big1/15tracks/figures/placefieldcorr_singleNeuronDirectly_orientation.pdf')
%%  fig 4e;population place field capacity
clear
savepath='/media/Big1/15tracks/placefieldCorr/';
load('/media/Big1/15tracks/placefieldCorr/PlFieldsCCMatrix.mat')

res=struct;n=0;
for irat=[1,2,3,4,5]
    seqID=out(irat).seqID;
    seqIDstr=out(irat).seqIDstr;
    
    for idir=1:2
        ind1=out(irat).ind_y{idir};
        seqID1=seqID(ind1,:);
        
        % same track, diff track
        inddiff=find( seqID1(:,strcmp( seqIDstr,'Embeded'))~=10 ...
            & seqID1(:,strcmp( seqIDstr,'Multitracks'))==0 ...
            & seqID1(:,strcmp( seqIDstr,'Repeat'))==1 ...
            & seqID1(:,strcmp(seqIDstr,'EntrDir'))==idir-1);
        indsame=find( seqID1(:,strcmp( seqIDstr,'Embeded'))~=10 ...
            & seqID1(:,strcmp( seqIDstr,'Repeat'))>1 ...
            & seqID1(:,strcmp(seqIDstr,'EntrDir'))==idir-1);
        
        indsamepair1=find( seqID1(:,strcmp( seqIDstr,'Track'))==1);
        indsamepair1=nchoosek(indsamepair1,2);
        indsamepair8=find( seqID1(:,strcmp( seqIDstr,'Track'))==8);
        indsamepair=[indsamepair1;indsamepair8(:)'];
        
        ndiff=length(inddiff);
        trackid=seqID1(inddiff,strcmp(seqIDstr,'Track'));
        ccr=out(irat).cc_y{idir};
        ccdiff=ccr(inddiff,inddiff);
        ccsame=[];
        
        for ic=1:size(indsamepair,1)
            ccsame=[ccsame;ccr(indsamepair(ic,1),indsamepair(ic,2))];
        end
        
        ccdiffmax=zeros(ndiff,1);
        for ic=1:ndiff
            %             ccdiffmax(ic)=max(ccdiff(ic,[1:ndiff]~=ic));
            tmp=ccdiff(1:ic,1:ic);tmp(tmp==1)=0;ccdiffmax(ic)=max(tmp(:));
        end
        
        ccdiffmaxRev=zeros(ndiff,1);
        ccdiffRev=ccdiff(end:-1:1,end:-1:1);
        for ic=1:ndiff
            tmp=ccdiffRev(1:ic,1:ic);tmp(tmp==1)=0;ccdiffmaxRev(ic)=max(tmp(:));
        end
        
        
        nrand=100;
        ccdiffmax_rand=zeros(ndiff,nrand);
        for irand=1:nrand
            randind=randperm(ndiff);
            ccdiff_rand=ccdiff(randind,randind);
            for ic=1:ndiff
                %             ccdiffmax(ic)=max(ccdiff(ic,[1:ndiff]~=ic));
                tmp=ccdiff_rand(1:ic,1:ic);tmp(tmp==1)=0;ccdiffmax_rand(ic,irand)=max(tmp(:));
            end
        end
        tmp=flipud(ccdiff); tmp(tmp==1)=NaN;
        diagm=diagmean(tmp);
        ccsamemean=mean(ccsame);
        
        % same maze, diff maze
        indmaze1=find( seqID1(:,strcmp( seqIDstr,'Embeded'))~=10 ...
            & seqID1(:,strcmp( seqIDstr,'Multitracks'))==0 ...
            & seqID1(:,strcmp( seqIDstr,'Repeat'))==1 ...
            & seqID1(:,strcmp( seqIDstr,'Track'))<8 ...
            & seqID1(:,strcmp(seqIDstr,'EntrDir'))==idir-1);
        indmaze2=find( seqID1(:,strcmp( seqIDstr,'Embeded'))~=10 ...
            & seqID1(:,strcmp( seqIDstr,'Multitracks'))==0 ...
            & seqID1(:,strcmp( seqIDstr,'Repeat'))==1 ...
            & seqID1(:,strcmp( seqIDstr,'Track'))>=8 ...
            & seqID1(:,strcmp(seqIDstr,'EntrDir'))==idir-1);
        ccmaze1=triu_kf(ccr(indmaze1,indmaze1));
        ccmaze2=triu_kf(ccr(indmaze2,indmaze2));
        ccdiffmaze=ccr(indmaze1,indmaze2);
        ccdiffmaze=ccdiffmaze(:);
        
        mse_maze=zeros(2,4);
        mse_maze_note='[mean,se],[maze1,maze2,diffmaze,samemaze]';
        mse_maze_Str={'maze1','maze2','diff','same'};
        [mse_maze(1,:),mse_maze(2,:)]=cellfun(@(x) mean_se(x,'mean'),{ccmaze1,ccmaze2,ccdiffmaze(:),[ccmaze1;ccmaze2]});
        pmaze12=ranksum(ccmaze1,ccmaze2);
        pmazesamediff=ranksum(ccdiffmaze(:),[ccmaze1;ccmaze2]);
        
        % multi track
        indmulti=find( seqID1(:,strcmp( seqIDstr,'Multitracks'))==1 ...
            & seqID1(:,strcmp(seqIDstr,'EntrDir'))==idir-1);
        ccmulti=ccr(indmulti(1),indmulti(2));
        
        %         res(irat,idir).ccdiffmax=ccdiffmax;
        %         res(irat,idir).ccsamemean=ccsamemean;
        %         res(irat,idir).ndiff=ndiff;
        %         res(irat,idir).trackid=trackid;
        %         res(irat,idir).
        EntrDir=idir-1;
        ratname=out(irat).ratname;
        excludeName={'out','res'};
        saveVar2structRes
        if isempty(fieldnames(res))
            res=res1;
        else
            res(irat,idir)=res1;
        end
    end
end
ccdiffmaxall=[cat(1,res(:,1).ccdiffmax),cat(1,res(:,2).ccdiffmax)];
ndiffx=cell2mat(cellfun(@(x) 1:x,{res(:,1).ndiff},'UniformOutput',false))';
% corr of cc and track#
trccRP=zeros(2,3);%[r,p][dir1,dir2,bothdir]
[trccRP(1,1),trccRP(2,1)]=corr(ccdiffmaxall(:,1),ndiffx);
[trccRP(1,2),trccRP(2,2)]=corr(ccdiffmaxall(:,2),ndiffx);
[trccRP(1,3),trccRP(2,3)]=corr(ccdiffmaxall(:),[ndiffx;ndiffx]);

mse_mazeall=cell(2,1);pmaze12=zeros(2,1);pmazesamediff=zeros(2,1);
for idir=1:2
    data={cat(1,res(:,idir).ccmaze1),cat(1,res(:,idir).ccmaze2),cat(1,res(:,idir).ccdiffmaze)...
        ,[cat(1,res(:,idir).ccmaze1);cat(1,res(:,idir).ccmaze2)]};
    [mse_mazeall{idir}(1,:),mse_mazeall{idir}(2,:)]=cellfun(@(x) mean_se(x,'median'),data);
    pmaze12(idir)=ranksum(data{1},data{2});
    pmazesamediff(idir)=ranksum(data{3},data{4});
end


%
f1=figure('position',[117          61        1606         572]);
for n=1:size(res,1)
    %     for idir=1:size(res,2)
    if isempty(res(n,1).ccdiffmax);continue;end
    subplot(1,size(res,1),n)
    idir=1;
    h1=plot(1:res(n,idir).ndiff,res(n,idir).ccdiffmax,'r-');hold on;
    h2=plot([1,res(n,idir).ndiff],res(n,idir).ccsamemean*[1,1],'rx-');
    h1r=plot(1:res(n,idir).ndiff,mean(res(n,idir).ccdiffmax_rand,2),'r:','color',h1.Color);
    idir=2;
    h3=plot(1:res(n,idir).ndiff,res(n,idir).ccdiffmax,'b-');hold on;
    h4=plot([1,res(n,idir).ndiff],res(n,idir).ccsamemean*[1,1],'bx-');
    h3r=plot(1:res(n,idir).ndiff,mean(res(n,idir).ccdiffmax_rand,2),'b:','color',h3.Color);
    set([h1;h2;h3;h4],'linewidth',2);
%     set(gca,'xtick',1:5:res(n,idir).ndiff,'xticklabel',num2strcell(res(n,idir).trackid(1:5:end)),...
%         'XTickLabelRotation',0);
    set(gca,'xtick',0:5:res(n,idir).ndiff);
    xlabel('track id');ylabel('Max CC ');
    title(sprintf('%s',res(n,1).ratname));
    %     end
end
setaxisformalAll(8)
legend([h1;h1r;h2;h3;h3r;h4],{'DiffTr Return','DiffTr Return Shuffle','SameTr Return',...
    'DiffTr EntrDir','DiffTr EntrDir Shuffle','SameTr EntrDir'},'location','best');

% select tracks
diagm=[]; tmp2=[];
for ii=1:length(res)
track_select=[1,2,3,4,8,9,10,11];
tmp=res(ii).ccdiff(track_select,track_select);
tmp(tmp==1)=NaN;
tmp2=cat(3,tmp2,tmp);
diagm=[diagm,diagmean(tmp,'LD')];
end
figure('position',[678    66   428   579]);
subplot(3,1,[1,2]);
imagesc(mean(tmp2,3));
set(gca,'xtick',1:length(track_select),'xticklabel',num2strcell(track_select),...
    'ytick',1:length(track_select),'yticklabel',num2strcell(track_select))
subplot(3,1,3);
plot(diagm);legend('location','eastoutside')
hold on;plot_errorbar(diagm');
diagm1=[repmat([1:size(diagm,1)]',5,1),diagm(:)];
% figure;plot(diagm1(:,1),diagm1(:,2),'k.');hold on; lsline
[dr,dp]=corr(diagm1(:,1),diagm1(:,2));
title(sprintf('r=%.3g, p=%.3g',dr,dp));
xlabel('Diagonal');
% 

f2=figure;
z=mean(cat(3,res(:).ccdiff),3);z(z==1)=NaN;
subplot(3,2,[1,3]);
imagesc(z);colorbar;title('pop vec corr');set(gca,'ydir','normal')
xlabel('track');ylabel('track');axis square;
setaxisformal
subplot(313);
z=cat(2,res(:).diagm);
zm=mean(z,2);
[zcc,zccp]=corr(zm,[1:length(zm)]');
plot(z);hold on;
plot(mean(z,2),'k','linewidth',2);
xlabel('time');ylabel('average along diagonal');set(gca,'xlim',[0,length(z)+1]);
title(sprintf('r=%.2g,p=%.2g',zcc,zccp));
setaxisformal(gca,12,0)
subplot(3,2,[2,4]);
plot(z,'color',0.5*[1,1,1]);hold on;
plot(mean(z,2),'k','linewidth',2);
xlabel('time');ylabel('average along diagonal');
title(sprintf('r=%.2g,p=%.2g',zcc,zccp));set(gca,'xlim',[0,length(z)+1]);
setaxisformal


f3=figure('position',[296    82   860   460]);
for idir=1:2
    subplot(1,3,idir);
    scatter(ndiffx,ccdiffmaxall(:,idir));
    lsline
    xlabel('#track');ylabel('Max CC ');
    title(sprintf('EntrDir%d,r=%.2g,p=%.2g',idir-1,trccRP(1,idir),trccRP(2,idir)));
end
subplot(1,3,3)
scatter([ndiffx;ndiffx],ccdiffmaxall(:));
lsline
xlabel('#track');ylabel('Max CC ');
title(sprintf('BothDir,r=%.2g,p=%.2g',trccRP(1,3),trccRP(2,3)));
setaxisformalAll(8)

% savepdf(f1,fullfile(savepath,'maxCCAllTr.pdf'));
% savepdf(f2,fullfile(savepath,'maxCCPositiveCorrAllTr.pdf'));
%  capacity estimation
f4=figure;
for idir=1:2
    ccdiff=cat(2,res(:,idir).ccdiffmax);
    ccdiffmax=mean(ccdiff,2);
    ccdiff_rand=cat(2,res(:,idir).ccdiffmax_rand);
    ccdiffmax_rand=mean(ccdiff_rand,2);
    
    ccsame=cat(2,res(:,idir).ccsame);
    ccsame=ccsame([1],:);
    [ccsamem,ccsamese]=mean_se(ccsame(:));
    ntrack=length(ccdiffmax);
    pval_sd=zeros(ntrack,1);
    for i=1:ntrack
        pval_sd(i)=ranksum(ccsame(:),ccdiff(i,:));
    end
    
    subplot(1,2,idir);
    
    plot(ccdiffmax,'o-');hold on;
    plot(ccdiffmax_rand,'x-');
    errorbar(1:length(ccdiffmax),ccsamem*ones(length(ccdiffmax),1),ccsamese*ones(length(ccdiffmax),1));
    %     sigstar1(1:ntrack,pval_sd,ccdiffmax);
    ylabel('max cc with previous track');xlabel('Track')
    title(['Track pair confusion (pop vec) Dir',num2str(idir)]);
    legend({'diff tracks','diff tracks shuffle time','same track'},'box','off','location','best');
end
setaxisformalAll

ccdiff=cat(2,res(:).ccdiffmax);
ccdiffmax=mean(ccdiff,2);
ccdiffRev1=cat(2,res(:).ccdiffmaxRev);
ccdiffRev=mean(cat(2,res(:).ccdiffmaxRev),2);
ccdiff_rand=cat(2,res(:).ccdiffmax_rand);
ccdiffmax_rand=mean(ccdiff_rand,2);

cc3={ccdiff(2:7,:),ccdiff_rand(2:7,:),ccdiffRev1(2:7,:)};
[cc3m,cc3se]=cellfun(@(x) mean_se(x(:)),cc3);
cc3p=[ranksum(cc3{1}(:),cc3{2}(:));ranksum(cc3{3}(:),cc3{2}(:))];

ccsame=cat(2,res(:).ccsame);
ccsame=ccsame([1],:);
[ccsamem,ccsamese]=mean_se(ccsame(:));
ntrack=length(ccdiffmax);
pval_sd=zeros(ntrack,1);
for i=1:ntrack
    pval_sd(i)=ranksum(ccsame(:),ccdiff(i,:));
end
f5=figure;
plot(ccdiffmax,'o-');hold on;
plot(ccdiffmax_rand,'x-');
errorbar(1:length(ccdiffmax),ccsamem*ones(length(ccdiffmax),1),ccsamese*ones(length(ccdiffmax),1));
%     sigstar1(1:ntrack,pval_sd,ccdiffmax);
ylabel('max cc with previous track');xlabel('Track')
title(['Track pair confusion (pop vec) Both Dir']);
legend({'diff tracks','diff tracks shuffle time','same track'},'box','off','location','best');
setaxisformalAll


f6=figure('color','w','position',[63         491        1463         820]);
subplot(221)
plot(ccdiffmax,'b-');hold on;
plot(ccdiffmax_rand,'b:');
% errorbar(1:length(ccdiffmax),ccsamem*ones(length(ccdiffmax),1),ccsamese*ones(length(ccdiffmax),1));
%     sigstar1(1:ntrack,pval_sd,ccdiffmax);
plot(ccdiffRev,'b-.');hold on;
ylabel('max cc with previous track');xlabel('Track')
title(['Track pair confusion (pop vec) Both Dir']);
legend({'diff tracks','diff tracks shuffle time','Rev order(tr15->tr1)'},'box','off','location','best');
subplot(222)
errorbarformal(1:3,cc3m,cc3se,{'diff tracks','diff tracks shuffle time','Rev order(tr15->tr1)'});
set(gca,'ylim',[0.1,0.4],'ytick',0:.1:1,'xticklabelrotation',30);
subplot(223)
plot(ccdiffmax-ccdiffmax_rand,'b-');hold on;
plot(ccdiffmax_rand-ccdiffmax_rand,'b:');
% errorbar(1:length(ccdiffmax),ccsamem*ones(length(ccdiffmax),1),ccsamese*ones(length(ccdiffmax),1));
%     sigstar1(1:ntrack,pval_sd,ccdiffmax);
plot(ccdiffRev-ccdiffmax_rand,'b-.');hold on;
ylabel('\DeltaConfusion');xlabel('Track')
title(['Diff']);
legend({'diff tracks','diff tracks shuffle time','Rev order(tr15->tr1)'},'box','off','location','best');

setaxisformalAll
%
ccdiffall=cat(3,res.ccdiff);
% parallel inner
ind1={[5,7];
[13,15
12,14]};
% perpendicular inner
ind2={[5,6; 6,7];
    [12,13;12,15;14,13;14,15]};

% ind3={[1,4;2,3];
%     [8,11;9,10]};
%  outer
ind_out={[1,4; 2,3],[8,11;9,10] ;
    [1,2;1,3;2,4;3,4],[8,9;8,10;9,11;10,11]};

ccout=cell(2,2); 
for i=1:2
    for m=1:2
    for j=1:size(ind_out{i},1)
        ccout{i,m}=[ccout{i,m};squeeze(ccdiffall(ind_out{i,m}(j,1),ind_out{i,m}(j,2),:))];
    end
    end
end
[ccoutm,ccoutse,ccoutp]=mean_se_cell(ccout,'mean');


cc57=cell(2,3);
for i=1:2
    for j=1:size(ind1{i},1)
        cc57{i,1}=[cc57{i,1};squeeze(ccdiffall(ind1{i}(j,1),ind1{i}(j,2),:))];
        
        cc57{i,2}=[cc57{i,2};ccdiff(ind1{i}(j,2),:)-ccdiff(ind1{i}(j,1),:)];
        cc57{i,3}=[cc57{i,3};ccdiffRev1(16-ind1{i}(j,1),:)-ccdiffRev1(16-ind1{i}(j,2),:)];
    end
end
[cc57m,cc57se,cc57p]=mean_se_cell(cc57(:,1));
[cc57maxm,cc57maxse,cc57maxp]=mean_se_cell(cc57(:,[2,3]));
f7=figure('position',[401   327   731   342]);
subplot(121);
errorbarformal(1:2,cc57m,cc57se,{'tr5-7','tr12-14,13-15'});
sigstar([1,2],cc57p(end))
ylabel('pop vector corr');
subplot(122);
errorbarformal(1:4,cc57maxm(:),cc57maxse(:),{'tr5-7','tr12-14,13-15','Revtr5-7','Revtr12-14,13-15'});
sigstar(mat2cell(cc57maxp(:,[1,2]),ones(size(cc57maxp,1),1),2),cc57maxp(:,end))
ylabel('max corr with previous track');
set(gca,'xticklabelrotation',30)
% eval_capacity_formula(ccdiffmax(1:end))
setaxisformalAll
 

cc56=cell(2,3);
for i=1:2
    for j=1:size(ind2{i},1)
        cc56{i,1}=[cc56{i,1};squeeze(ccdiffall(ind2{i}(j,1),ind2{i}(j,2),:))];
        
        cc56{i,2}=[cc56{i,2};ccdiff(ind2{i}(j,2),:)-ccdiff(ind2{i}(j,1),:)];
        cc56{i,3}=[cc56{i,3};ccdiffRev1(16-ind2{i}(j,1),:)-ccdiffRev1(16-ind2{i}(j,2),:)];
    end
end
[cc56m,cc56se,cc56p]=mean_se_cell(cc56(:,1));
[cc56maxm,cc56maxse,cc56maxp]=mean_se_cell(cc56(:,[2,3]));

load(fullfile('/media/Big1/15tracks/placefieldCorr/','PlFieldsCCMatrixAlldir_geometryFeatures.mat'),'datacc'...
    );
paramaze12={datacc{4}{4},datacc{5}{4}}; % paralell maze1; paralell maze2
[pallm,pallse]=mean_se_cell(paramaze12');
pallp=ranksum(paramaze12{1},paramaze12{2});

%
% parallel inner
ind13={[5,7]
[12,14]
[13,15]
};
% perpendicular inner
ind23={[5,6];
    [6,7];
    [12,13];
    [13,14;
    12,15]
    [14,15]};

cc53para=cell(length(ind13),1);
    for j=1:length(ind13)
        cc53para{j}=[cc53para{j};squeeze(ccdiffall(ind13{j}(1),ind13{j}(2),:))];
    end
[cc53m,cc53se,cc53p]=mean_se_cell(cc53para);
[cc53tmpy,cc53tmpx]=cell2group(cc53para);
[cc53r,cc53rp]=corr(cc53tmpx,cc53tmpy);
%  test_multigroup_pairwise(cell2mat(cc53para'))
cc54orth=cell(length(ind23),1);
    for j=1:length(ind23)
        for jj=1:size(ind23{j},1)
        cc54orth{j}=[cc54orth{j};squeeze(ccdiffall(ind23{j}(jj,1),ind23{j}(jj,2),:))];
        end
    end
[cc54m,cc54se,cc54p]=mean_se_cell(cc54orth);
[cc54tmpy,cc54tmpx]=cell2group(cc54orth);
[cc54r,cc54rp]=corr(cc54tmpx,cc54tmpy);


% -------- each individual track pair
f534=figure('position',[401   327   731   342]);
subplot(121);
plot_errorbar_with_data(cc53tmpy,cc53tmpx,{'5-7','12-14','13-15'});
sigstar(mat2cell(cc53p(:,[1,2]),ones(3,1),2),cc53p(:,end))
ylabel('pop vector corr');xlabel('Track pair');
title(sprintf('r=%.4g, p=%.4g',cc53r,cc53rp));
subplot(122);
plot_errorbar_with_data(cc54tmpy,cc54tmpx,{'5-6','6-7','12-13','12-14','13-14','14-15'});
sigstar(mat2cell(cc54p(:,[1,2]),ones(size(cc54p,1),1),2),cc54p(:,end))
ylabel('pop vector corr');xlabel('Track pair');
setaxisformalAll
title(sprintf('r=%.4g, p=%.4g',cc54r,cc54rp));


f8=figure('position',[654    94   638   427]);
data1=[paramaze12(:);cc56(:,1);cc57(:,1)];
hb=errorbarformal(1:6,[pallm;cc56m;cc57m],[pallse;cc56se;cc57se],...
    {'Maze1 paral.','Maze2 paral.','Maze1 Inner perp.','Maze2 Inner perp.','Maze1 Inner paral.','Maze2 Inner paral.'});
set(hb(1),'CData',[1,0,0; 0,0,1;1,0,0; 0,0,1;1,0,0; 0,0,1;]);
set(hb(1),'FaceColor','Flat');
for i1=1:length(data1)
    plot(i1+(rand(size(data1{i1}))-0.5)/8,abs(data1{i1}),'k.');
end

sigstar({[1,2],[3,4],[5,6]},[pallp,cc56p(end),cc57p(end)]);hold on;

% errorbarformal(3:4,cc57m,cc57se,{'tr5-7','tr12-14,13-15'});
% sigstar([3,4],cc57p(end))
ylabel('pop vector corr');
set(gca,'xticklabelrotation',30);
title(sprintf('p=%.2g, %.2g, %.2g',pallp,cc56p(end),cc57p(end)));
% eval_capacity_formula(ccdiffmax(1:end))
setaxisformalAll

f9=figure('position',[654    94   638   427]);
data1=[ccout(2,:)';ccout(1,:)';cc56(:,1);cc57(:,1)];
hb=errorbarformal(1:8,[ccoutm(2,:)';ccoutm(1,:)';cc56m;cc57m],[ccoutse(2,:)';ccoutse(2,:)';cc56se;cc57se],...
    {'Maze1 Out perp.','Maze2 Out perp.','Maze1 Out paral.','Maze2 Out paral.','Maze1 Inner perp.','Maze2 Inner perp.','Maze1 Inner paral.','Maze2 Inner paral.'});
set(hb(1),'CData',[1,0,0; 0,0,1;1,0,0; 0,0,1;1,0,0; 0,0,1;1,0,0; 0,0,1;]);
set(hb(1),'FaceColor','Flat');
for i1=1:length(data1)
    plot(i1+(rand(size(data1{i1}))-0.5)/8,abs(data1{i1}),'k.');
end

sigstar({[1,2],[3,4],[5,6],[7,8]},[ccoutp(6,end),ccoutp(1,end),cc56p(end),cc57p(end)]);hold on;

% errorbarformal(3:4,cc57m,cc57se,{'tr5-7','tr12-14,13-15'});
% sigstar([3,4],cc57p(end))
ylabel('pop vector corr');
set(gca,'xticklabelrotation',30);
title(sprintf('p=%.3g, %.3g, %.3g, %.3g',ccoutp(6,end),ccoutp(1,end),cc56p(end),cc57p(end)));
% eval_capacity_formula(ccdiffmax(1:end))
setaxisformalAll

% savepdf(f1,'/media/Big1/15tracks/figures/pop_vec__maxcc_eachrat.pdf')
% savepdf(f2,'/media/Big1/15tracks/figures/pop_vec_maxcc_timeDynamic.pdf')
% savepdf(f6,'/media/Big1/15tracks/figures/pop_vec_maxcc_cumulative_timeDynamic.pdf')
% savepdf(f7,'/media/Big1/15tracks/figures/pop_vec_maxcc_paraleltrack_time.pdf')
%
% savepdf(f8,'/media/Big1/15tracks/figures/pop_vec_maxcc_perpendicular_track_time.pdf')
%  savepdf(f9,'/media/Big1/15tracks/figures/pop_vec_maxcc_innerouter_track_time.pdf')
 
%% fig s4d. half session rate map corr
load('/media/Big1/15tracks/RunSeqCorr/1stHalf2ndHalf_rateCorr.mat');
%
f1=figure;
subplot(211);
histogram(rp(:,1),0:.1:1,'normalization','probability')
xlabel('Rate map correlation');ylabel('Frequency');
subplot(212);
for irat=1:length(ccrp)
    plot(ccrp(irat).rateccrp(:,1),'-','color',0.5*[1,1,1]);
    hold on;
end
plot(1:15,rpm,'k','linewidth',2);
xlabel('Track');ylabel('Rate map correlation');
f2=figure;
subplot(1,3,[1,2]);
errorshade21(1:15,rpm,rpse);
xlabel('Track');ylabel('Rate map correlation');
set(gca,'ylim',[0,1],'ytick',0:.5:1);
subplot(1,3,3);
boxplot(rp(:,1),'Notch','off','Symbol','');
% boxplot(rp(:,1),rati(:,1),'Notch','off','Symbol','');
hold on;plot([0,2],[0,0],'r:');
sigstar1(1,rpp,1);
ylabel('Rate map correlation');set(gca,'ylim',[0,1],'ytick',0:.5:1);
setaxisformalAll(12,0);

[m1,se1]=mean_se(rp(:,1))

figure;
histogram(rp_tr1(:,1),0:.1:1);
ylabel('Rate map correlation of track1rep1rep2');

% savepdf(f2,'/media/Big1/15tracks/figures/halfsession_rateCorr.pdf')
%% fig 4b 4c; laps within each track discrminability
clear
 load('/media/Big1/15tracks/placefieldCorr/EachlapDiscriminationAll.mat','out','snrAll');
% load('/media/Big1/15tracks/placefieldCorr/EachlapDiscriminationAll_occupancy.mat','out','snrAll');

%
seqIDstr=out(1).seqIDstr;
seqID1=cat(1,snrAll.seqID);
ccrp=cat(1,snrAll.ccrp);
ccr0p=signrank(ccrp(:,1),0);
ccmse=zeros(1,2);[ccmse(1),ccmse(2)]=mean_se(ccrp(~isnan(ccrp(:,1)),1));
ccedge=-1:.2:1;ccedgex=edge2x(ccedge);
ccrpN=histcounts(ccrp(:,1),ccedge);
ccrpNsig=histcounts(ccrp(ccrp(:,2)<.05,1),ccedge);
[ccrpm,ccrpse]=mean_se(remove_nan_vector( ccrp(:,1)));
signrank(ccrp(:,1))


[sum(ccrp(:,1)<0&ccrp(:,2)<.05),sum(ccrp(:,1)<0&ccrp(:,2)>.05),...
    sum(ccrp(:,1)>0&ccrp(:,2)>.05),sum(ccrp(:,1)>0&ccrp(:,2)<.05)]
ccrpsigind=find(ccrp(:,2)<.05);
ccrppositiveind=find(ccrp(:,1)>0);


ccrpRand=cat(1,snrAll.ccrpShuf); 
ccrpRand=reshape(permute(ccrpRand,[1,3,2]),size(ccrpRand,1)*size(ccrpRand,3),size(ccrpRand,2));
ccrpNrand=histcounts(ccrpRand(:,1),ccedge); ccrpNrand=ccrpNrand/sum(ccrpNrand);
ccrpN0=ccrpN/sum(ccrpN);

% figure;histogram(ccrpRand(:,1))
[~,p]=kstest2(ccrp(:,1),ccrpRand(:,1))
p=ranksum(ccrp(:,1),ccrpRand(:,1))
% signrank(ccrp(:,1))

f12=figure('color','w');
ylimdata=[0,max(ccrpN0)*1.75];
hb=bar(ccedgex,[ccrpNrand],1,'facecolor',0.5*[1,1,1],'linestyle','none');
hold on;
hb1=bar(ccedgex,[ccrpN0],1,'facecolor','b','linestyle','none');
set([hb,hb1],'FaceAlpha',0.5);
h1=plot(ccedgex,ccrpN0,'color','b');
h2=plot(ccedgex,ccrpNrand,'color',0.5*[1,1,1]);
plot([0,0],ylimdata,'k:');
xlabel('Correlation(lap~discrim)');ylabel('Frequency');
legend([hb1,hb],{'Data','Lap shuffle'},'location','best');
ylim(ylimdata)
setaxisformal

% figure;for i=1:size(seqID1,2);subplot(3,4,i);histogram(seqID1(ccrpsigind,i));title(seqIDstr{i});end
% suptitle('Significant positive tracks');
% figure;for i=1:size(seqID1,2);subplot(3,4,i);histogram(seqID1(ccrppositiveind,i));title(seqIDstr{i});end
% suptitle('Positive tracks');

f1=figure('color','w','position',[  513    96   695   394]);
ylimdata=[0,max(ccrpN)*1.1];
subplot(121);
h1=bar(ccedgex,ccrpN,.5,'w');hold on;
h2=bar(ccedgex,ccrpNsig,.5,'b');
plot([0,0],ylimdata,'k:');
xlabel('Correlation(lap~discrim)');ylabel('# tracks');
legend({'All','Signicficant'},'location','best');
ylim(ylimdata)
setaxisformal
subplot(122);
errorbarformal(1,ccmse(1),ccmse(2));
hold on;
plot([0,2],[0,0]);
sigstar1(1,ccr0p,ccmse(1));

% distributionPlot(ccrp(:,1))

ylabel('Correlation');
setaxisformal




% example discrimination

for i=7%1:length(out)
    
    tx=cell(out(i).ntrack,1);
    for itr=1:out(i).ntrack
        s=cell2mat(out(i).spktime(:,itr));
        lapx=unique(s(:,6));nlap=length(lapx);
        tx{itr}=zeros(nlap,1);
        for ilap=1:nlap
            tx{itr}(ilap)=min(s(s(:,6)==lapx(ilap),1));
        end
        %         tx{itr}( tx{itr}==0)=[];
    end
    f2=figure('color','w');
    nrow=ceil(sqrt(out(i).ntrack));
    for itr=9%1:out(i).ntrack
        %     subplot(nrow,nrow,itr);
        y=snrAll(i).snr{itr}; if isempty(y);continue;end
        plot(1:length(y),y,'.');
        xlim([0,length(y)+1]);ylim([min(y),max(y)*1.1])
        title(sprintf('%dR%.1gP%.1g',out(i).trackx(itr),...
            snrAll(i).ccrp(itr,1),snrAll(i).ccrp(itr,2)));
        xlabel('lap');ylabel('Discriminability');
    end
    setaxisformalAll
    % savepdf(gcf,fullfile('/media/Big1/15tracks/placefieldCorr',...
    %     [out(i).ratname,'EntrDir',num2str(out(i).entrdir),'DiscrimCorrWithlaps.pdf']));
    % close all
end


% savepdf(f12,'/media/Big1/15tracks/figures/lap_discrim_distribution2bar.pdf')
% 
% savepdf(f1,'/media/Big1/15tracks/figures/lap_discrim_distribution.pdf')
% savepdf(f2,'/media/Big1/15tracks/figures/lap_discrim_example.pdf')


% savepdf(gcf,fullfile('/media/Big1/15tracks/placefieldCorr',...
%     ['DiscrimCorrWithlapsCCsummary.pdf']));
% save('/media/Big1/15tracks/placefieldCorr/EachlapDiscriminationAll.mat','out','snrAll');
%

%% fig s4a. example demo for discriminability
clear;close all;
load('/media/Big1/15tracks/placefieldCorr/EachlapDiscriminationAll_occupancy_2dir.mat','out');
iout=1


spktime=out(iout).spktime; %1,x;2y; 3linear distance; 4 speed;5dir;6dir angular;7time
spkhist=out(iout).spkhist;
tracksession=out(iout).tracksession;
times=out(iout).times;
trackdir=out(iout).trackdir;
dirValue=out(iout).dirValue;
pf=out(iout).pf;
shortlap_ind=out(iout).shortlap_ind;

iseq=1;
for ilap=[9,22]

cellind=[5,6,26,23];
spktime1=spktime(cellind,iseq);
vlim=5;
load(fullfile(out(iout).matpath,'Maze.mat'),'times','mazel') % lap time

posvel=mazel{tracksession(iseq)}(:,[1 2 4 7 8 5 6]); 
%1,x;2y; 3linear distance; 4 speed;5dir;6dir angular;7time      
time1lap=times{tracksession(iseq),trackdir(iseq)};

goodlap=find(shortlap_ind{iseq}==0);

                pos1=posvel( posvel(:,5)==dirValue(trackdir(iseq)) ...
                    & posvel(:,6)~=0 & posvel(:,7)>time1lap(goodlap(ilap),1) ...
                    & posvel(:,7)<time1lap(goodlap(ilap),2),:);
               
linecolor=jet(length(cellind));

figure;
subplot(3,2,1)
plot(pos1(:,3),ones(size(pos1(:,7))),'k','linewidth',2);hold on;
for ic=1:length(cellind)
   thislap=spktime{cellind(ic)}(spktime{cellind(ic)}(:,3)>vlim ...
       & spktime{cellind(ic)}(:,6)==goodlap(ilap),:); 
   % 1spike time, 2position (linear), 3velocity (cm/s),
        %       4direction/lap 5direction (from angular velocity)
   plot(thislap(:,2),ones(size(thislap(:,1))),'.','markersize',12,'color',linecolor(ic,:));   
end
title(ilap)
subplot(3,2,[3,5]);
for ic=1:length(cellind)
    thislap=squeeze(spkhist{iseq}(ilap,cellind(ic),:)); thislap=thislap/max(thislap);
    plot(1:length(thislap),thislap,'color',linecolor(ic,:)); hold on;
end
xlabel('location');ylabel('neuron');title('lap activity');
subplot(1,5,[4,5]);
for iseq1=1:15
    for ic=1:length(cellind)
    thistr=pf{iseq1}(cellind(ic),:); thistr=thistr/max(thistr);
    plot(linspace(0,1,length(thistr)),thistr+iseq1,'color',linecolor(ic,:)); hold on;
    end
end
xlabel('location');ylabel('track');title('track activity');
end
saveallfig2pdf('/media/Big1/15tracks/figures/discriminability_Demo.pdf');

% figure;plot_plfield(pf{iseq});
%%  discriminability early late
clear
load('/media/Big1/15tracks/placefieldCorr/EachlapDiscrimination_earlylate.mat');

snr4=cellfun(@(x) mean(x),cat(1,snrAll.snr));

f1=figure;
% subplot(121);
[m,se]=mean_se(snr4);
[out4,groupPair,pval]=test_multigroup_pairwise(snr4);
errorbarformal(1:2,m,se,snr_str);
sigstar(groupPair,pval);
ylabel('discriminability');
title(sprintf('p=%.2g',pval));
setaxisformal;

% savepdf(f1,'/media/Big1/15tracks/figures/discriminability_earlyLate.pdf')
%% fig s4b. Orientation selectivity early late
clear
load('/media/Big1/15tracks/placefieldCorr/EachlapDiscrimination_earlylate.mat');


osi1=cat(2,osiAll.osi)';
goodind=cat(2,osiAll.goodind)';
goodind=goodind & (osi1(:,1)>0.4&osi1(:,2)>0.4);
osi1good=osi1(goodind,:);
osi1bad=osi1(~goodind,:);

figure;
subplot(131);
[~,p]=plot_scatter(osi1);
title(sprintf('All neuron \n p=%.2g',p))
xlabel('OSI(early)');ylabel('OSI(late)');
subplot(132);
[~,p]=plot_scatter(osi1good);
title(sprintf('ori selective \n p=%.2g',p))
xlabel('OSI(early)');ylabel('OSI(late)');
subplot(133);
[~,p]=plot_scatter(osi1bad);
title(sprintf('ori non-select\n p=%.2g',p))
xlabel('OSI(early)');ylabel('OSI(late)');

f2=figure('position',[  109   191   937   285]);
osi_str={'early','late'};
subplot(131);
[m,se]=mean_se(osi1);
p=signrank(osi1(:,1),osi1(:,2));
[ydata,gr]=cell2group({osi1(:,1),osi1(:,2)}');
% errorbarformal(1:2,m,se);
plot_errorbar_with_data(ydata,gr);
hold on;
sigstar([1,2],p);
title(sprintf('All neuron \n p=%.2g',p))
set(gca,'xtick',1:2,'xticklabel',osi_str,'xticklabelrotation',30);
ylabel('OSI');
subplot(132);
[m,se]=mean_se(osi1good);
p=signrank(osi1good(:,1),osi1good(:,2));
% errorbarformal(1:2,m,se);hold on;
[ydata,gr]=cell2group({osi1good(:,1),osi1good(:,2)}');
plot_errorbar_with_data(ydata,gr);
hold on;
sigstar([1,2],p);
title(sprintf('ori selective \n p=%.2g',p))
set(gca,'xtick',1:2,'xticklabel',osi_str,'xticklabelrotation',30,'ylim',[0.4,1]);
ylabel('OSI');
subplot(133);
[m,se]=mean_se(osi1bad);
p=signrank(osi1bad(:,1),osi1bad(:,2));
% errorbarformal(1:2,m,se);hold on;
[ydata,gr]=cell2group({osi1bad(:,1),osi1bad(:,2)}');
plot_errorbar_with_data(ydata,gr);
hold on;
sigstar([1,2],p);
title(sprintf('ori non-selective \n p=%.2g',p))
set(gca,'xtick',1:2,'xticklabel',osi_str,'xticklabelrotation',30);
ylabel('OSI');
setaxisformalAll

% savepdf(f2,'/media/Big1/15tracks/figures/OSI_earlyLate.pdf')
savepdf(f2,'/media/Big1/15tracks/figures/OSI_earlyLate_dot.pdf')

%% fig s4c; run speed
% edit runspeed_each_session_15track.m
rawdatapath='/media/Big1/15tracks/';
ratnames={'jsrat1','jsrat2','jsrat3','jsrat4','jsrat5'};
    
savepath=['/media/Big1/15tracks/runSpeed'];
% mkdir(savepath);

v=cell(length(ratnames),1);
verr=cell(length(ratnames),1);
for irat=1:length(ratnames)
    cd(fullfile(rawdatapath,ratnames{irat}));
    irat
    load('Maze.mat','mazel');
% 1 mean x, 2 mean y, 3 angle, 4 distance in cm!,
%  5 direction sign, 6 timestamp, 7 velocity
    [v{irat},verr{irat}]=cellfun(@(x) mean_se(x(x(:,7)~=0,7),'median'),mazel);
end

v=cat(2,v{:});
verr=cat(2,verr{:});
vm=mean(v');
%
nses=12;
figure;
for i=1:length(ratnames)
    subplot(2,3,i);
shadedErrorBar(1:nses,v(:,i),verr(:,i),'b');
title(ratnames{i});
xlabel('run session');ylabel('speed(cm/s)');
set(gca,'ylim',[0,40]);
end

figure;
plot(v,'color',0.5*[1,1,1]);hold on;
plot(vm,'b','linewidth',2);
xlabel('Run session');ylabel('Run speed (cm/s)');
set(gca,'ylim',[0,30],'xlim',[1,length(vm)]);
setaxisformal

savepdf(gcf,'/media/Big1/15tracks/figures/runspeed_across_session.pdf');
%%  discriminability early late orientation

clear
load('/media/Big1/15tracks/placefieldCorr/EachlapDiscrimination_earlylate.mat');
snr4=cellfun(@(x) mean(x),cat(1,snrAll.snr));
snr4All=snr4;

% f1=figure;
% subplot(121);
% [m,se]=mean_se(snr4);
% [out4,groupPair,pval]=test_multigroup_pairwise(snr4);

load('/media/Big1/15tracks/placefieldCorr/EachlapDiscrimination_earlylateOrientation.mat');

snr6=cat(2,snr4All,snr4);
gr=repmat(1:size(snr6,2),size(snr6,1),1);
gr=gr(:);
ydata=snr6(:);

f1=figure('position',[154   508   922   420]);
subplot(121);
% [m,se]=mean_se(snr4);
% [out4,groupPair,pval]=test_multigroup_pairwise(snr4);
% errorbarformal(1:4,m,se,snr_str);
% sigstar(groupPair,pval);
% ylabel('discriminability');
[hdot,hbar,ydatam,ydatase]=plot_errorbar_with_data(ydata,gr,{'early','late','earlySame','lateSame','earlyDiff','lateDiff'});
[out6,groupPair,pval]=test_multigroup_pairwise(snr6);
sigstar(groupPair,pval);
set(gca,'ylim',[0,20]);
ylabel('discriminability');

subplot(122);
[m1,se1]=mean_se(snr4d);
[out4,groupPair,pval1]=test_multigroup_pairwise(snr4d);
errorbarformal(1:4,m1,se1,{'sameori increment','diffori increment','diff-same early','diff-same late'});
sigstar(groupPair,pval1);
ylabel('discriminability ratio');
setaxisformalAll
disp(m1)
table(groupPair,pval1)


snr4sameDdiff=[snr4(:,1)./snr4(:,3),snr4(:,2)./snr4(:,4)];
pval_samediff=signrank(snr4sameDdiff(:,1),snr4sameDdiff(:,2));
figure;
[m,se]=mean_se(snr4sameDdiff);
errorbarformal(1:2,m,se,{'early same/diff','late same/diff'});hold on;
sigstar([1,2],pval_samediff);
ylabel('discriminability ratio');
%

% savepdf(f1,'/media/Big1/15tracks/figures/discrim_orientation_earlyLate.pdf')
% savepdf(f1,'/media/Big1/15tracks/figures/discrim_orientation_earlyLate_dot.pdf')
%% miss assignment to previous or future tracks
clear
load('/media/Big1/15tracks/placefieldCorr/EachlapDiscriminationAll_occupancy_2dir.mat','out');

correctRate15=cat(3,out.correctRate15);

missAssign=[];c=[];
for i=1:length(out)
    y=[sum(tril_kf(correctRate15(:,:,i))),sum(triu_kf(correctRate15(:,:,i)))];
    missAssign=[missAssign;y/15];
    
    c=[c;diag(correctRate15(:,:,i))];
    
end


figure;
subplot(121);
histogram(c,0.7:0.01:1);
% plot_cdf_with_meanse({c},{'Correct'})
xlabel('Correct rate(%)');
subplot(122)
plot(missAssign'*100,'o-');
set(gca,'xtick',[1,2],'xticklabel',{'previous','future'},'xlim',[0,3],'xticklabelrotation',30);
ylabel('Proportion of laps (%)');
 setaxisformalAll
 
% savepdf(gcf,'/media/Big1/15tracks/figures/missAssign_to_previous_or_future_tracks.pdf')
%% fig.4a; laps within each track discrminability--orientation
clear
load('/media/Big1/15tracks/placefieldCorr/EachlapDiscriminationAllOrientation.mat');

f1=figure;
plot_cdf_with_meanse(snrsamediff,{'sameori','diffori'});
set(gca,'xscale','log');title(sprintf('p=%.2g',snrsamediff_pval));
xlabel('Discriminability');ylabel('Frequency');
setaxisformal

% savepdf(f1,'/media/Big1/15tracks/figures/lap_discrim_orientation.pdf')
%% number of neuron in each rat
pathall={
    '/media/Big1/15tracks/jsrat1/';
    '/media/Big1/15tracks/jsrat2/';
    '/media/Big1/15tracks/jsrat3/';
    '/media/Big1/15tracks/jsrat4/';
    '/media/Big1/15tracks/jsrat5/';
    };
ncell=zeros(length(pathall),2);
for ipath=1:length(pathall)
    matpath=pathall{ipath};
    ratname=ratnameFromMatpath(matpath);
    % pyramidal cell
    [Cl1, Cl2] =textread(fullfile(matpath,'celltype'), '%s %s');
    Pyr=strcmp(Cl1,'p');
    ncell(ipath,1)=sum(Pyr);
    ncell(ipath,2)=sum(Pyr==0);
end
[m1,se1]=mean_se(ncell(:,1));
[m2,se2]=mean_se(ncell(:,2));
f1=figure;
histogram(ncell(:,1),0:5:80);hold on;
histogram(ncell(:,2),0:5:80);
legend({'Pyramidal','Interneuron'});
setaxisformalAll;
xlabel('Number of neurons');ylabel('Number of animals');
title(sprintf('Pyr %.2g+1%.2g,Interneu %.2g+1%.2g',m1,se1,m2,se2));
% savepdf(f1,['/media/Big1/15tracks/figures/Numberofneuron_eachrat.pdf'])
